Correlation workshop

Micom 2021

This is an informal and hopefully intuitive way to teach the nature of correlations. Please consider textbooks and review-articles for detailed study. This project can be downloaded from GitHub or using git directly using shell commands:

git clone https://github.com/bioinformatics-leibniz-hki/correlation-workshop

We can deploy a Docker container containing RStudioServer with all required tools using GNU make:

cd correlation-workshop
make

The prompt will display the URL in which the environment can be accessed. The username is rstudio with password password.

Computer Science doctorates vs. arcade gaming revenue

Correlation is way to quantify the relationship between two variables. This relationship is monotonous and often even linear. We start with looking at a dataset from the U.S. Census Bureau and National Science Foundation showing the number of computer science doctorated and the revenue of arcade games in billion USD (“Spurious Correlations,” n.d.).

First, we load the data using functions from the R tidyverse packages:

library(tidyverse)

arcades <- read_csv("data/arcades.csv")

arcades %>% summary()

      year        doctorates        revenue     
 Min.   :2000   Min.   : 809.0   Min.   :1.176  
 1st Qu.:2002   1st Qu.: 862.5   1st Qu.:1.247  
 Median :2004   Median :1038.5   Median :1.371  
 Mean   :2004   Mean   :1195.1   Mean   :1.442  
 3rd Qu.:2007   3rd Qu.:1571.5   3rd Qu.:1.641  
 Max.   :2009   Max.   :1787.0   Max.   :1.803  

Let’s plot one variable against the other and do a linear regression:

arcades %>%
  ggplot(aes(doctorates, revenue)) +
  stat_smooth(method = "lm") +
  geom_point() +
  labs(
    x = "Computer science doctorates",
    y = "Arcade revenue (billion USD)"
  )

Linear correlation

The higher the number of computer science doctorates, the higher is the revenue in the arcade gaming market! Correlation values are always ranging between -1 and 1. The value is negative if one variable continuously increases whereas the other decreases and positive if both variables are either increasing or decreasing. The value is zero if there is no monotonous relationship.

Let’s calculate the linear Pearson’s product-moment correlation:

cor.test(~ doctorates + revenue, data = arcades, method = "pearson")

    Pearson's product-moment correlation

data:  doctorates and revenue
t = 16.182, df = 8, p-value = 2.138e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.935915 0.996586
sample estimates:
      cor 
0.9850653 

The correlation coefficient is estimated to be \(0.985\) indicating a strong positive relationship which fits with the plot. Most basic statistical tests, including those for correlations, are actually nothing else than linear models (Lindeløv, n.d.). The estimated slope of the fitted curve is the strength of the correlation. We have to normalize our data using z-scaling (\(z(x)=\frac{x-\mu}{\sigma}\)) to ensure that the estimated slope is independent of the unit and lies between -1 and 1. This will give us values with a mean of 0 and a variance of 1.

arcades_scaled <-
  arcades %>%
  mutate(
    doctorates = doctorates %>% scale(),
    revenue = revenue %>% scale()
  )

arcades_scaled %>%
  pivot_longer(c(doctorates, revenue)) %>%
  group_by(name) %>%
  summarise(mean = mean(value), variance = var(value))

# A tibble: 2 x 3
  name            mean variance[,1]
  <chr>          <dbl>        <dbl>
1 doctorates  2.41e-16           1.
2 revenue    -1.85e-16           1 

Now we can do a linear regression on one variable using the other as a covariate to get the same correlation coefficient. We ignore the intercept here because this tell us nothing about the relationship between the two variables and should be almost 0 due to the scaling:

lm(formula = doctorates ~ 1 + revenue, data = arcades_scaled) %>% summary()

Call:
lm(formula = doctorates ~ 1 + revenue, data = arcades_scaled)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.27327 -0.12554 -0.00275  0.12641  0.29887 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 5.565e-16  5.775e-02    0.00        1    
revenue     9.851e-01  6.088e-02   16.18 2.14e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1826 on 8 degrees of freedom
Multiple R-squared:  0.9704,    Adjusted R-squared:  0.9666 
F-statistic: 261.8 on 1 and 8 DF,  p-value: 2.138e-07

Ranked correlation

Often, we do not require the relationship to be linear and proportional. For instance, this is the case if we want to get the correlation of the abundance of a species with time in an exponential growing setting. Now we don’t care about the size of the gradient triangles anymore. We can remove this effect by ranking the data and doing the linear model again:

arcades_ranked <-
  arcades %>%
  mutate(
    doctorates = doctorates %>% rank(),
    revenue = revenue %>% rank()
  )


arcades_ranked %>%
  ggplot(aes(doctorates, revenue)) +
  stat_smooth(method = "lm") +
  geom_point() +
  coord_fixed() +
  labs(
    x = "Computer science doctorates (rank)",
    y = "Arcade revenue (rank)"
  )

lm(formula = doctorates ~ 1 + revenue, data = arcades_ranked) %>% summary()

Call:
lm(formula = doctorates ~ 1 + revenue, data = arcades_ranked)

Residuals:
     Min       1Q   Median       3Q      Max 
-3.12727 -0.02121  0.25455  0.68182  1.21212 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.4667     0.8843   0.528 0.612017    
revenue       0.9152     0.1425   6.421 0.000204 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.295 on 8 degrees of freedom
Multiple R-squared:  0.8375,    Adjusted R-squared:  0.8172 
F-statistic: 41.23 on 1 and 8 DF,  p-value: 0.0002045

The test is called Spearman’s rank correlation:

cor.test(~ doctorates + revenue, data = arcades, method = "spearman")

    Spearman's rank correlation rho

data:  doctorates and revenue
S = 14, p-value = 0.0004667
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.9151515 

Note that the linear correlation coefficient is higher indicating that this relationship is indeed linear: \(r_{pearson} = 0.985 > r_{spearman} = 0.912\)

Partial correlation

One might ask at this point why there is a relationship between the number of computer science doctorates and the arcade revenue. This question is totally valid and the example dataset was chosen on purpose: To illustrate the ubiquitous phenomena called spurious correlation. Pearson himself already noted this potential for false positive discoveries (Pearson 1897). This is especially the case if a third confounding variable is influencing both variables in the same way. Here, this is obviously just a coincidence so the third confounding variable is the time itself.

To mitigate this issue, we control for this third variable by fitting both other variables against the year. The residuals of the linear model is what can’t be explained by time so we can do the correlation test on the residuals:

doctorates_residuals <- lm(doctorates ~ year, data = arcades_scaled) %>% residuals()
revenue_residuals <- lm(revenue ~ year, data = arcades_scaled) %>% residuals()

cor.test(~ doctorates_residuals + revenue_residuals, method = "pearson")

    Pearson's product-moment correlation

data:  doctorates_residuals and revenue_residuals
t = 6.2546, df = 8, p-value = 0.0002445
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6604151 0.9790925
sample estimates:
      cor 
0.9111654 

The linear Pearson correlation is now indeed lower indicating that time is a confounding variable and the relationship is just a coincidence. This process of doing the test only on a fraction of the total variance (namely the residuals) is called partial correlation. There is also the function ppcor::pcor.test doing this:

library(ppcor)
pcor.test(
  x = arcades_scaled$doctorates,
  y = arcades_scaled$revenue,
  z = arcades_scaled$year,
  method = "pearson"
)

   estimate      p.value statistic  n gp  Method
1 0.9111654 0.0006301684  5.850676 10  1 pearson

Local correlation

There are relationships between two variables which are statistical depended but can not be discovered using correlation test. Correlation implies statistical dependency but not vice versa. This is especially the case in which the two variables are not monotonously linked but in a cyclic way. Let’s simulate the abundance of some predator and some prey species to estimate the relationship of hunting:

hunting <-
  tibble(time = seq(from = 0, to = 3, length.out = 50)) %>%
  mutate(
    predator = sin(time),
    prey = 1 + cos(time - 0.5)
  )

hunting %>%
  pivot_longer(
    cols = c(predator, prey),
    names_to = "animal",
    values_to = "abundance"
  ) %>%
  ggplot(aes(time, abundance, color = animal)) +
  geom_point()

The pearson correlation is almost zero here:

cor.test(~ predator + prey, data = hunting, method = "pearson")

    Pearson's product-moment correlation

data:  predator and prey
t = 0.68411, df = 48, p-value = 0.4972
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1851467  0.3665860
sample estimates:
       cor 
0.09826513 

This is because the relationship between the predator and the prey is not linear:

hunting %>%
  ggplot(aes(predator, prey)) +
  geom_point()

But we can still do a correlation test if we just look at a small time window in which the relationship between the predator and the prey can be assumed to be linear:

hunting_local <- hunting %>% filter(time > 1 & time < 2)

hunting_local %>%
  pivot_longer(
    cols = c(predator, prey),
    names_to = "animal",
    values_to = "abundance"
  ) %>%
  ggplot(aes(time, abundance, color = animal)) +
  geom_point() +
  stat_smooth(method = "lm")

cor.test(~ predator + prey, data = hunting_local, method = "pearson")

    Pearson's product-moment correlation

data:  predator and prey
t = -1.7183, df = 14, p-value = 0.1078
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.75653210  0.09881308
sample estimates:
       cor 
-0.4173342 

Now we can observe the negative local correlation between the predator and the prey. An elaborated version of this idea is called Local similarity analysis (LSA) (Ruan et al. 2006)

Species and pathway abundance in healthy and disesed cohorts

Let’s continue with a more realistic dataset. We have healthy subjects (Nielsen et al. 2014) and patients suffering from Colorectal cancer (CRC) (Zeller et al. 2014). This dataset contains whole metagenomic sequencing data which was processed by Humann2 to get abundance profiles for bacterial species and pathways. Data was retrieved using the R-package curatedMetagenomicData (Pasolli et al. 2017).

features <- read_csv("data/features.csv")
samples <- read_csv("data/samples.csv")
lineages <- read_csv("data/lineages.csv")

samples %>% mutate_if(is.character, as.factor) %>% summary()

   sample_id               sample_alias           project   body_site 
 S01    : 1   CCIS12656533ST-4-0 : 1    NielsenHB_2014:16   stool:40  
 S02    : 1   CCIS16383318ST-4-0 : 1    ZellerG_2014  :24             
 S03    : 1   CCIS29688262ST-20-0: 1                                  
 S04    : 1   CCIS41288781ST-4-0 : 1                                  
 S05    : 1   CCIS48725289ST-4-0 : 1                                  
 S06    : 1   CCIS53355328ST-4-0 : 1                                  
 (Other):34   (Other)            :34                                  
    disease        age        age_category    gender   country 
 CRC    :20   Min.   :29.00   adult :29    female:10   DEU:11  
 healthy:20   1st Qu.:50.75   senior:11    male  :14   DNK:11  
              Median :59.00                NA's  :16   ESP: 5  
              Mean   :58.98                            FRA:13  
              3rd Qu.:67.25                                    
              Max.   :90.00                                    
                                                               

The R-package ComplexHeatmap can be used to create a heatmap showing the abundance profile. Similar samples and species with a high Spearman correlation are displayed in neighboring rows.

library(ComplexHeatmap)

features_mat <-
  features %>%
  pivot_longer(-sample_id, names_to = "feature", values_to = "abundance") %>%
  left_join(lineages) %>%
  
  group_by(feature) %>%
  filter(kingdom == "Bacteria" & var(abundance) > 0) %>%
  
  group_by(sample_id) %>%
  mutate(abundance = abundance / sum(abundance)) %>%

  # long table to wide matrix
  ungroup() %>%
  select(sample_id, feature, abundance) %>%
  pivot_wider(names_from = feature, values_from = abundance) %>%
  select(-sample_id) %>%
  as.matrix()

Heatmap(
  matrix = features_mat,
  name = "Abundance (TSS)",

  column_title = "Species",
  clustering_distance_columns = "spearman",
  show_column_names = FALSE,

  row_title = "Sample",
  clustering_distance_rows = "spearman",

  right_annotation = HeatmapAnnotation(
    which = "row",
    Disease = samples$disease,
    Age = samples$age,
    Country = samples$country,
    Gender = samples$gender
  )
)

Sparse correlation

Count data in ecology has often many zeros. This zero-inflation violates the assumption of a linear model. It assumes a normal distribution of the residuals (Gauss–Markov theorem). The majority of the taxa are not present in an average sample and only in a few samples.

Methods like SparCC take this into account (Friedman and Alm 2012). We use the function correlate_spiec_easi_sparcc inside the source directory which is based on the implementation in SpiecEasi (Layeghifard, Hwang, and Guttman 2018). This allows us to get an object of class tidygraph with one table for the edges and nodes.

# load  libraries and functions
source("src/setup.R")

features_sparcc <-
  lineages %>%
  filter(kingdom == "Bacteria") %>%
  pull(feature) %>%
  head(20)

mat_sparcc <-
  features %>%
  pivot_longer(-sample_id, names_to = "feature", values_to = "abundance") %>%
  left_join(lineages) %>%
  filter(feature %in% features_sparcc) %>%
  # long table to wide matrix
  select(sample_id, feature, abundance) %>%
  pivot_wider(names_from = feature, values_from = abundance, values_fill = list(abundance = 0)) %>%
  select(-sample_id) %>%
  as.matrix()

res_sparcc <- correlate_spiec_easi_sparcc(
  data = mat_sparcc,
  nodes = lineages,
  iterations = 10,
  bootstraps = 20
  )

res_sparcc

# A tbl_graph: 11 nodes and 8 edges
#
# An unrooted forest with 3 trees
#
# Node Data: 11 x 13 (active)
  feature feature_type kingdom phylum order class family genus species degree
  <chr>   <chr>        <chr>   <chr>  <chr> <chr> <chr>  <chr> <chr>    <dbl>
1 Eubact… taxon        Bacter… Firmi… Clos… Clos… Eubac… Euba… Eubact…      1
2 Bacter… taxon        Bacter… Bacte… Bact… Bact… Bacte… Bact… Bacter…      1
3 Alisti… taxon        Bacter… Bacte… Bact… Bact… Riken… Alis… Alisti…      1
4 Eubact… taxon        Bacter… Firmi… Clos… Clos… Eubac… Euba… Eubact…      2
5 Bifido… taxon        Bacter… Actin… Acti… Bifi… Bifid… Bifi… Bifido…      1
6 Rumino… taxon        Bacter… Firmi… Clos… Clos… Lachn… Blau… Rumino…      1
# … with 5 more rows, and 3 more variables: component <int>, closeness <dbl>,
#   betweeness <dbl>
#
# Edge Data: 8 x 31
   from    to estimate p.value q.value from_feature from_feature_ty…
  <int> <int>    <dbl>   <dbl>   <dbl> <chr>        <chr>           
1     3     8    0.329       0       0 Alistipes p… taxon           
2     8    10   -0.159       0       0 Bacteroides… taxon           
3     8     9    0.370       0       0 Bacteroides… taxon           
# … with 5 more rows, and 24 more variables: from_kingdom <chr>,
#   from_phylum <chr>, from_order <chr>, from_class <chr>, from_family <chr>,
#   from_genus <chr>, from_species <chr>, from_degree <dbl>,
#   from_component <int>, from_closeness <dbl>, from_betweeness <dbl>,
#   to_feature <chr>, to_feature_type <chr>, to_kingdom <chr>, to_phylum <chr>,
#   to_order <chr>, to_class <chr>, to_family <chr>, to_genus <chr>,
#   to_species <chr>, to_degree <dbl>, to_component <int>, to_closeness <dbl>,
#   to_betweeness <dbl>

Let’s plot the network and use functions of the package ggraph and correlate_plot from this repository:

library(ggraph)

res_sparcc %>%
  correlate_plot() +
  geom_node_point(aes(color = phylum), size = 10) +
  geom_node_text(aes(label = feature))

Compositional correlation

Sequencing data is compositional: Instead of counting the species in a absolute way, one just have proportions of the reads assigned to a particular taxon. This can have severe implications (Gloor et al. 2017). Imagine one growing species and many constant ones in a longitudinal dataset. The relative abundance of all other taxa are decreasing, because the library size (total number of reads sequenced) is fixed. This happens for all other taxa in the same way introducing lots of false positive correlations. BAnOCC is a tool which takes this compositionality into account (Schwager et al. 2017). The sum for all abundances must be 1 for each sample. Let’s use their example dataset to get the correlations:

library(banocc)
data(compositions_null)

correlate_banocc(compositions_null)

Running /usr/local/lib/R/bin/R CMD SHLIB foo.c
gcc -I"/usr/local/lib/R/include" -DNDEBUG   -I"/usr/local/lib/R/site-library/Rcpp/include/"  -I"/usr/local/lib/R/site-library/RcppEigen/include/"  -I"/usr/local/lib/R/site-library/RcppEigen/include/unsupported"  -I"/usr/local/lib/R/site-library/BH/include" -I"/usr/local/lib/R/site-library/StanHeaders/include/src/"  -I"/usr/local/lib/R/site-library/StanHeaders/include/"  -I"/usr/local/lib/R/site-library/RcppParallel/include/"  -I"/usr/local/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/usr/local/lib/R/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fpic  -g -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -g  -c foo.c -o foo.o
In file included from /usr/local/lib/R/site-library/RcppEigen/include/Eigen/Core:88,
                 from /usr/local/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13,
                 from <command-line>:
/usr/local/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:1: error: unknown type name ‘namespace’
  613 | namespace Eigen {
      | ^~~~~~~~~
/usr/local/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:613:17: error: expected ‘=’, ‘,’, ‘;’, ‘asm’ or ‘__attribute__’ before ‘{’ token
  613 | namespace Eigen {
      |                 ^
In file included from /usr/local/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
                 from /usr/local/lib/R/site-library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13,
                 from <command-line>:
/usr/local/lib/R/site-library/RcppEigen/include/Eigen/Core:96:10: fatal error: complex: No such file or directory
   96 | #include <complex>
      |          ^~~~~~~~~
compilation terminated.
make: *** [/usr/local/lib/R/etc/Makeconf:167: foo.o] Error 1

# A tbl_graph: 0 nodes and 0 edges
#
# An empty graph
#
# Node Data: 0 x 4 (active)
# … with 4 variables: degree <dbl>, component <int>, closeness <dbl>,
#   betweeness <dbl>
#
# Edge Data: 0 x 12
# … with 12 variables: from <int>, to <int>, estimate <dbl>, conf_alpha <dbl>,
#   from_degree <dbl>, from_component <int>, from_closeness <dbl>,
#   from_betweeness <dbl>, to_degree <dbl>, to_component <int>,
#   to_closeness <dbl>, to_betweeness <dbl>

A comprehensive titorial is also available.

Same correlation

This dinosaur dataset tells us that one can create almost any pattern having the same statistical values including the correlation (Tran 2020; Autodesk 2017):

Which correlation?

Correlation detection strategies in microbial data sets vary widely in sensitivity and precision. An overview is given in (Weiss et al. 2016).

Drake workflow

Analyses usually consists if multiple steps which are depended on each other. We can write this work flow as a tree and use the work flow management R package drake. This has several advantages:

• running multiple steps (called targets in drake) in parallel
• Saving intermediate results
• Rerun only targets which actually changed

correlation_plan <- drake_plan(
  features = file_in("data/features.csv") %>% read_csv(),
  lineages = file_in("data/lineages.csv") %>% read_csv,
  feature_types = c("taxon", "pathway"),
  cor_methods = c("spearman", "pearson"),
  
  matricies = target({
    matrix <-
      features %>%
      pivot_longer(-sample_id, names_to = "feature", values_to = "abundance") %>%
      left_join(lineages) %>%
      
      # get current feature type
      filter(feature_type == feature_types) %>%
      select(sample_id, feature, abundance) %>%
      
      # subsetting
      sample_frac(0.05) %>%
      
      # long table to wide matrix
      pivot_wider(
        names_from = feature,
        values_from = abundance,
        values_fill = list(abundance = 0)
      ) %>%
      select(-sample_id) %>%
      as.matrix()
    
    tibble(feature_type = feature_types, matrix = list(matrix))
  },
  dynamic = map(feature_types)
  ),
  
  correlations = target({
    correlation <- correlate(data = matricies$matrix[[1]], method = cor_methods)
    
    tibble(
      feature_type = feature_types,
      cor_method = cor_methods,
      correlation = list(correlation)
    )
  },
  dynamic = cross(feature_types, cor_methods)
  )
)

correlation_plan %>% make()

Results were stored on hard disk. They can be loaded into the current R session:

loadd(correlations)

correlations

# A tibble: 4 x 3
  feature_type cor_method correlation
  <chr>        <chr>      <list>     
1 taxon        spearman   <tbl_grph> 
2 taxon        pearson    <tbl_grph> 
3 pathway      spearman   <tbl_grph> 
4 pathway      pearson    <tbl_grph> 

A more elaborated work flow is part of this repository (See report). It includes also the part of getting the data. Drake plans are tibbles of steps and can be plotted:

plan %>% plot()

The targets can be made by running the script analyze.R. It is recommended to put this as a background job of RStudio to run it in an isolated environment.

rstudioapi::jobRunScript("analyze.R")

Ensemble correlation

Often, an edge is only considered if it is significant in multiple correlation methods. These ensemble approaches tend to be more robust having less false positive interactions:

loadd(coabundances)

graphs <-
  coabundances %>%
  filter(feature_type == "pathway" & disease == "CRC") %>%
  select(method, coabundance) %>%
  deframe()

ensemble_coabundance(graphs$spearman, graphs$banocc, method = "intersection")

# A tbl_graph: 17 nodes and 18 edges
#
# An undirected simple graph with 4 components
#
# Node Data: 17 x 5 (active)
  feature                                  degree component closeness betweeness
  <chr>                                     <dbl>     <int>     <dbl>      <dbl>
1 adenosine ribonucleotides de novo biosy…      2         3   0.00417          0
2 adenine and adenosine salvage III             1         1   0.00538          0
3 5-aminoimidazole ribonucleotide biosynt…      2         3   0.00417          0
4 6-hydroxymethyl-dihydropterin diphospha…      1         1   0.00543          0
5 5-aminoimidazole ribonucleotide biosynt…      2         3   0.00417          0
6 4-deoxy-L-threo-hex-4-enopyranuronate d…      4         1   0.00562          8
# … with 11 more rows
#
# Edge Data: 18 x 14
   from    to estimate_rcorr multiedges from_feature from_degree from_component
  <int> <int>          <dbl>      <int> <chr>              <dbl>          <int>
1     6    10          0.8            2 4-deoxy-L-t…           4              1
2    10    12          0.753          2 (5Z)-dodec-…           3              1
3     8    17          0.711          2 acetylene d…           4              2
# … with 15 more rows, and 7 more variables: from_closeness <dbl>,
#   from_betweeness <dbl>, to_feature <chr>, to_degree <dbl>,
#   to_component <int>, to_closeness <dbl>, to_betweeness <dbl>

No correlation but covariance

Correlation assumes that the pair of two variables are conditionally independend from all others. But often, an interaction is only there in a presence or absence if a third variable. SpiecEasi tries to account for this (Kurtz et al. 2015). An application of this is finding relationships between features of different types. Let’s load both species and pathways from the Colorectal Cancer (CRC) cohort from the main plan in this repository (20 samples and 20 features each).

loadd(inter_feature_type_mats)

matricies_crc <-
  inter_feature_type_mats %>%
  filter(disease == "CRC") %>%
  pull(mat) %>%
  first()

matricies_crc %>% map(dim)

[[1]]
[1] 20 20

[[2]]
[1] 20 20

Let’s estimate covariance relationships:

crc_spiec_easi <- matricies_crc %>% correlate_spiec_easi(method = "mb", nodes = lineages)
crc_spiec_easi

# A tbl_graph: 27 nodes and 32 edges
#
# An undirected multigraph with 8 components
#
# Node Data: 27 x 13 (active)
  feature feature_type kingdom phylum order class family genus species degree
  <chr>   <chr>        <chr>   <chr>  <chr> <chr> <chr>  <chr> <chr>    <dbl>
1 Faecal… taxon        Bacter… Firmi… Clos… Clos… Rumin… Faec… Faecal…      4
2 Paraba… taxon        Bacter… Bacte… Bact… Bact… Porph… Para… Paraba…      2
3 Alisti… taxon        Bacter… Bacte… Bact… Bact… Riken… Alis… Alisti…      2
4 Dorea … taxon        Bacter… Firmi… Clos… Clos… Lachn… Dorea Dorea …      2
5 Rumino… taxon        Bacter… Firmi… Clos… Clos… Lachn… Blau… Rumino…      1
6 Clostr… taxon        Bacter… Firmi… Clos… Clos… Clost… Clos… Clostr…      2
# … with 21 more rows, and 3 more variables: component <int>, closeness <dbl>,
#   betweeness <dbl>
#
# Edge Data: 32 x 29
   from    to estimate from_feature from_feature_ty… from_kingdom from_phylum
  <int> <int>    <dbl> <chr>        <chr>            <chr>        <chr>      
1    18    19   0.300  (5Z)-dodec-… pathway          <NA>         <NA>       
2    25    26   0.0877 acetyl-CoA … pathway          <NA>         <NA>       
3    10    26  -0.131  adenine and… pathway          <NA>         <NA>       
# … with 29 more rows, and 22 more variables: from_order <chr>,
#   from_class <chr>, from_family <chr>, from_genus <chr>, from_species <chr>,
#   from_degree <dbl>, from_component <int>, from_closeness <dbl>,
#   from_betweeness <dbl>, to_feature <chr>, to_feature_type <chr>,
#   to_kingdom <chr>, to_phylum <chr>, to_order <chr>, to_class <chr>,
#   to_family <chr>, to_genus <chr>, to_species <chr>, to_degree <dbl>,
#   to_component <int>, to_closeness <dbl>, to_betweeness <dbl>

crc_spiec_easi %>% correlate_plot()

These are the relationships found between a species and a pathway:

crc_spiec_easi %>%
  activate(edges) %>%
  as_tibble() %>%
  filter(from_feature_type != to_feature_type) %>%
  select(from_feature, to_feature, estimate)

# A tibble: 5 x 3
  from_feature               to_feature                                 estimate
  <chr>                      <chr>                                         <dbl>
1 Faecalibacterium prausnit… 4-deoxy-L-threo-hex-4-enopyranuronate deg…   0.0347
2 Ruminococcus obeum         acetyl-CoA fermentation to butanoate II      0.0477
3 Alistipes shahii           adenosine nucleotides degradation II        -0.0947
4 Alistipes shahii           adenosine nucleotides degradation II        -0.0461
5 Faecalibacterium prausnit… 4-deoxy-L-threo-hex-4-enopyranuronate deg…   0.0354

Network analysis

These correlation networks were created by the main drake work flow:

loadd(coabundances)

coabundances

# A tibble: 14 x 4
   disease method   feature_type coabundance
   <chr>   <chr>    <chr>        <list>     
 1 CRC     mb       all          <tbl_grph> 
 2 healthy mb       all          <tbl_grph> 
 3 healthy spearman taxon        <tbl_grph> 
 4 healthy sparcc   taxon        <tbl_grph> 
 5 healthy banocc   taxon        <tbl_grph> 
 6 CRC     spearman taxon        <tbl_grph> 
 7 CRC     sparcc   taxon        <tbl_grph> 
 8 CRC     banocc   taxon        <tbl_grph> 
 9 healthy spearman pathway      <tbl_grph> 
10 healthy sparcc   pathway      <tbl_grph> 
11 healthy banocc   pathway      <tbl_grph> 
12 CRC     spearman pathway      <tbl_grph> 
13 CRC     sparcc   pathway      <tbl_grph> 
14 CRC     banocc   pathway      <tbl_grph> 

The tidygraph package provides functions to calculate node and edge attributes based on the topological structure:

coabundances$coabundance[[1]] %>%
  activate(nodes) %>%
  mutate(
    # number of edges connected to a node
    degree = centrality_degree(),
    
    # Kleinberg's hub centrality scores
    hub = centrality_hub()
  ) %>%
  as_tibble() %>%
  select(feature, degree, hub)

# A tibble: 9 x 3
  feature                                         degree   hub
  <chr>                                            <dbl> <dbl>
1 Alistipes shahii                                     1 0.187
2 adenosine ribonucleotides de novo biosynthesis       1 0.375
3 adenine and adenosine salvage III                    2 0    
4 5-aminoimidazole ribonucleotide biosynthesis II      3 0.938
5 5-aminoimidazole ribonucleotide biosynthesis I       3 1.00 
6 adenosine nucleotides degradation II                 2 0.500
7 (5Z)-dodec-5-enoate biosynthesis                     2 0    
8 8-amino-7-oxononanoate biosynthesis I                2 0    
9 4-aminobutanoate degradation V                       2 0    

Objects created by correlate are shipped with basic topology metrics. We can extract them to make plot for comparison:

coabundances %>%
  filter(feature_type == "all") %>%
  select(-feature_type) %>%
  mutate(nodes = coabundance %>% map(~ .x %>% activate(nodes) %>% as_tibble)) %>%
  select(-coabundance) %>%
  unnest(nodes) %>%
  pivot_longer(cols = c(degree, closeness, betweeness), names_to = "metric") %>%

  ggplot(aes(feature_type, value)) +
      geom_boxplot() +
      facet_wrap(~metric, scales = "free_y") +
      labs(x = "Feature type", y = "Topology score")

Or compare taxonomic correlations between healthy and diseased cohorts:

coabundances %>%
  filter(feature_type == "taxon" & method == "sparcc") %>%
  select(-feature_type) %>%
  mutate(nodes = coabundance %>% map(~ .x %>% activate(nodes) %>% as_tibble)) %>%
  select(-coabundance) %>%
  unnest(nodes) %>%
  pivot_longer(cols = c(degree, closeness, betweeness), names_to = "metric") %>%
  
  ggplot(aes(disease, value)) +
    geom_boxplot() +
    ggpubr::stat_compare_means() +
    facet_wrap(~metric, scales = "free_y")

The graph with all attributes of edges and nodes can be saved as a file e.g. to visualize it in other tools like Cytoscape:

library(igraph)

coabundances %>%
  filter(
    disease == "healthy" &
      feature_type == "pathway" &
      method == "spearman"
  ) %>%
  pull(coabundance) %>%
  first() %>%
  as.igraph() %>%
  write.graph("healthy_pathway_spearman.graphml", format = "graphml")

References

Autodesk. 2017. “Same Stats, Different Graphs: Generating Datasets with Varied Appearance and Identical Statistics Through Simulated Annealing.” https://www.autodesk.com/research/publications/same-stats-different-graphs.

Friedman, Jonathan, and Eric J Alm. 2012. “Inferring Correlation Networks from Genomic Survey Data.” PLoS Comput Biol 8 (9): e1002687.

Gloor, Gregory B, Jean M Macklaim, Vera Pawlowsky-Glahn, and Juan J Egozcue. 2017. “Microbiome Datasets Are Compositional: And This Is Not Optional.” Frontiers in Microbiology 8: 2224.

Kurtz, Zachary D, Christian L Müller, Emily R Miraldi, Dan R Littman, Martin J Blaser, and Richard A Bonneau. 2015. “Sparse and Compositionally Robust Inference of Microbial Ecological Networks.” PLoS Comput Biol 11 (5): e1004226.

Layeghifard, Mehdi, David M Hwang, and David S Guttman. 2018. “Constructing and Analyzing Microbiome Networks in R.” In Microbiome Analysis, 243–66. Springer.

Lindeløv, Jonas Kristoffer. n.d. “Common Statistical Tests Are Linear Models (or: How to Teach Stats).” https://lindeloev.github.io/tests-as-linear/.

Nielsen, H Bjørn, Mathieu Almeida, Agnieszka Sierakowska Juncker, Simon Rasmussen, Junhua Li, Shinichi Sunagawa, Damian R Plichta, et al. 2014. “Identification and Assembly of Genomes and Genetic Elements in Complex Metagenomic Samples Without Using Reference Genomes.” Nature Biotechnology 32 (8): 822–28.

Pasolli, E., L. Schiffer, P. Manghi, A. Renson, V. Obenchain, D. T. Truong, F. Beghini, et al. 2017. “Accessible, curated metagenomic data through ExperimentHub.” Nat Methods 14 (11): 1023–4.

Pearson, Karl. 1897. “Mathematical Contributions to the Theory of Evolution.—on a Form of Spurious Correlation Which May Arise When Indices Are Used in the Measurement of Organs.” Proceedings of the Royal Society of London 60 (359-367): 489–98.

Ruan, Quansong, Debojyoti Dutta, Michael S Schwalbach, Joshua A Steele, Jed A Fuhrman, and Fengzhu Sun. 2006. “Local Similarity Analysis Reveals Unique Associations Among Marine Bacterioplankton Species and Environmental Factors.” Bioinformatics 22 (20): 2532–8.

Schwager, Emma, Himel Mallick, Steffen Ventz, and Curtis Huttenhower. 2017. “A Bayesian Method for Detecting Pairwise Associations in Compositional Data.” PLoS Computational Biology 13 (11): e1005852.

“Spurious Correlations.” n.d. http://www.tylervigen.com/spurious-correlations.

Tran, Khuyen. 2020. “Can Datasets of a Dinosaur and a Circle Have Identical Statistics?” https://towardsdatascience.com/how-to-turn-a-dinosaur-dataset-into-a-circle-dataset-with-the-same-statistics-64136c2e2ca0.

Weiss, Sophie, Will Van Treuren, Catherine Lozupone, Karoline Faust, Jonathan Friedman, Ye Deng, Li Charlie Xia, et al. 2016. “Correlation Detection Strategies in Microbial Data Sets Vary Widely in Sensitivity and Precision.” The ISME Journal 10 (7): 1669–81.

Zeller, Georg, Julien Tap, Anita Y Voigt, Shinichi Sunagawa, Jens Roat Kultima, Paul I Costea, Aurélien Amiot, et al. 2014. “Potential of Fecal Microbiota for Early-Stage Detection of Colorectal Cancer.” Molecular Systems Biology 10 (11): 766.